Appendix C: Elevation Data
Overview
The goal of this tutorial is to introduce the packages and techniques that are commonly utilized to manipulate raster data and create grid cells. We will use the elevation data of several larges counties in the Midwest as an example, and special emphasis will be placed on introducing the functionalities of the “velox” package. In short, our objectives are to:
- Learn the features of the “velox” package
- Inspect the elevation data of several large counties in the Midwest
- Output a new shapefile that contains the elevation data for each grid
Environment Setup
Input/Output
Our inputs include a shapefile of several large counties in the Midwest, a grid data file, and a raster file that contains the elevation data which is the focus of this tutorial. The files can be found here. The files are available for download on this project’s GitHub repo.
Our output is a new shapefile in which we write the the elevation data for each grid.
Load Libraries
We will use the following packages in this tutorial:
sf: to conduct basic spatial data manipulationtidyverse: to perform simple statistical analysesraster: to manipulate and analyze raster datagstat: to conduct geostatistical modeling and simulationtamp: to create spatial data visualizationvelox: to manipulate raster data in time efficient manner
The only package that is worth special mention the “velox,” which has excellent performance in fast raster data manipulation. Unfortunately, this package is no longer available. For the lack of substitutes, however, we will still utilize the velox package in this tutorial.
To install velox from archive, we can run the following code.
library(devtools)
devtools::install_github("https://github.com/hunzikp/velox")(Note, you may have to install an installer, such as Xcode, in order to compile the package locally. Some helpful information can be found here.)
library(sf)
library(tidyverse)
library(raster)
library(gstat)
library(tmap)
library(velox)Load Data
Loading data is the next step. We load the “grid” data (in kilometers), which is named km.grid, the shapefile of large counties, which is named lac, and also the elevation data, which is named lac.elevation.
km.grid <- st_read("./data/Km_Grid")## Reading layer `Km_Grid' from data source `/Users/LorenzMenendez/Desktop/OpenAirQ-toolkit/data/Km_Grid' using driver `ESRI Shapefile'
## Simple feature collection with 29903 features and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -88.78167 ymin: 40.73704 xmax: -86.46887 ymax: 43.19439
## geographic CRS: WGS 84lac <- st_read("./data/LargeAreaCounties")## Reading layer `LargeAreaCounties' from data source `/Users/LorenzMenendez/Desktop/OpenAirQ-toolkit/data/LargeAreaCounties' using driver `ESRI Shapefile'
## Simple feature collection with 21 features and 8 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -88.77879 ymin: 40.73641 xmax: -86.46629 ymax: 43.21831
## geographic CRS: WGS 84lac.elevation <- raster("./data/lac.elevation.grd")The elevation data comes from the 3D Elevation Program (3DEP) of United States Geological Survey. One nice thing about this program is that all of its data are available “free of charge and without use restrictions”. If you are interested in learning more about what data are available, feel free to explore this website.
For more detailed description of the data, please refer to the main chapters of this tutorial book.
Inspecting Data
Before we proceed any further, let’s check the resolution of our elevation data first. The resolution is shown in degrees because that is the unit of CRS.
It is worth pointing out that the resolution of the original raster data file is finer than what is shown below. Unfortunately, because Github only supports data files that are less than 100 MB, we downsampled the raster data to decrease its size. In other words, we decreased the resolution of the raster data by applying the aggregate function in the “raster” package. The good news is that the current resolution (after the downsampling) is sufficient for our purposes.
raster::res(lac.elevation)## [1] 0.0005555556 0.0005555556We can now plot a raster map, as shown below.
tm_shape(lac.elevation) +
tm_raster(alpha = .5) +
tm_shape(lac) +
tm_borders() +
tm_layout(legend.outside = TRUE,
main.title = "Elevation of Large Counties in Midwest",
main.title.size = 1.2)Data Manipulation
Raster data manipulation is where the velox packages really shines as it contains a large number of functions that you will find useful in handling raster data. For example, aggregate allows users to aggregate a VeloxRaster object to a lower resolution. Another function, crop, as its name suggests, lets you crop a VeloxRaster function. If you want to rasterize a set of polygons, then you can use rasterize. There are many more functions that are of interest in raster data manipulation, but it is impossible to list all of them here. Detailed documentation about this package can be found here, and you should feel free to explore it at your own pace. In our tutorial, we will mainly be using velox and extract.
First, we call the velox function to convert lac.elevation into a VeloxRaster object. Then we extract the mean value of each grid (km \(\times\) km), and we name it km.elecation.vx. Note that “velox” is indeed good at fast raster data manipulation - running the two lines of code takes very little time.
elevation.vx <- velox(lac.elevation)
km.elevation.vx <- elevation.vx$extract(km.grid, fun = mean)head(km.elevation.vx)## [,1]
## 1 266.9327
## 2 267.1034
## 3 267.6558
## 4 264.9626
## 5 272.6876
## 6 285.0676Let’s pause for a second and examine the object km.elevation.vx. As we would expect, it is just a list of the averages of each grid cell.
Next, we join this elevation data with km.grid, which we loaded earlier but haven’t touched at all.
km.grid$Elevation <- as.numeric(km.elevation.vx) Plotting Grid Map
With the data ready, we can now plot a graph.
tm_shape(km.grid) +
tm_fill("Elevation") +
tm_layout(legend.outside = TRUE, main.title = "Elevation of Large Counties in Midwest", main.title.size = 1.2)